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ABSTRACT 


Most first-order upstream conservative differencing methods can capture 
shocks quite well for one -dimensional problems. A direct application of 
these first-order methods to two-dimensional problems does not necessarily 
produce the same type of accuracy unless the shocks are locally aligned with 
the mesh. Harten has recently developed a second-order high-resolution ex- 
plicit method for the numerical computation of weak solutions of 
one-dimensional hyperbolic conservation laws. The main objectives of this 
paper are (a) to examine the shock resolution of Harten's method for a 
two-dimensional shock reflection problem, (b) to study the use of a 
high-resolution scheme as a post-processor to an approximate steady-state so- 
lution, and (c) to construct an implicit method in the delta-form using Har- 
ten's scheme for the explicit operator and a simplified iteration matrix for 



the implicit operator 
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I. INTRODUCTION 


The stability, accuracy and efficiency of a numerical scheme should have 
a symbiotic relationship with the types of equations that one Intends to 
solve numerically and the methods for treating shocks, contact discontinui- 
ties, and numerical boundary conditions. A typical complaint about the ap- 
plication of a scheme that is developed under the guideline of linear theory 
for the compressible Euler or Navier-Stokes equations is that the resulting 
scheme is not robust and or not accurate enough. Often the root of the ins- 
tability and inaccuracy lies in the improper treatment of nonlinear effects 
[1,2] (shocks and contact discontinuities) and numerical boundary condition 
procedures [3|. 

In problems with shocks where conventional second or higher-order accu- 
rate central spatial difference methods are used, the resulting numerical so- 
lution exhibits overshoots and undershoots in the vicinity of discontinuities 
[4]. The oscillations not only degrade the accuracy but can cause nonlinear 
instabilities. One remedy is to add numerical dissipation. Unfortunately, 
ad hoc methods of adding dissipation generally smear the discontinuities. 

This work was motivated by Marten's [5] recent success in developing a 
high-resolution second-order explicit method for one-dimensional hyperbolic 
conservation laws which has the following properties: 

(ci) -’me is developed in conservation form to ensure that the 


limit is a weak solution. 
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(b) The scheme satisfies a proper entropy Inequality [1-2] to ensure 
that the limit solution will have only physically-relevant discontinuities. 

(c) The scheme Is designed such that the amount of numerical dissipa- 
tion used produces highly accurate weak solutions. 

The goal of this paper Is to examine the shock resolution of Harten's 
method for a two-dimensional gas-dynamic problem, and to investigate other 
£ replications of his method including the possible extension to a 
high-resolution implicit method for both one- and two-dimension problems. 

We have applied Harten's method to shock tube [5] and 
quasi-one-dimensional nozzle problems. Also, we have made a direct applica- 
tion of his method to two-dimensional transient [5] and steady-state 
gas-dynamic problems. In both one and two dimensions, accurate weak solu- 
tions were obtained. 



A brief description of a particular version of Harten's method and its 
extension is given in sections II and III. In section IV, we will describe 
the two-dimensional implementation of Harten's method by a fractional step 
method and a modified implicit method. Numerical experiments for quasi-one 
dimensional nozzle problems ard a two-dimensional shock reflection problem 
are given in sections V and VI. 


II. HARTEN'S SCHEME 
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To aid the development and motivation of the next section, we briefly 
describe a particular version of Marten's second-order explicit scheme for a 
one-dlmenslonal system of conservation laws. The reader should refer to the 
original paper [5] for a more detailed description. 

Consider a system of conservation laws of the form 


9t dx 


(1) 



where U is the vector of m conserved variables and F is the flux vector. The 
Jacobian matrix A(U) = ?F(U)/ 5U has real eigenvalues (V" > , • . . , 

and a complete set of (right) eigenvectors. Let R =» (R^ ,R^, . . . ,R"*) be a ma- 
trix whose columns are the right eigenvectors of A and let L be a matrix 
whose rows are the left eigenvectors of A normalized such t it LR ■ I. 


Let IJ. = V(U ■ denote an average of U- and U-., . A simple exam- 

J.+ 72 f 

pie is the arithmetic average )/2 (see Roe [6] for a more 

sophisticated formula for for inviscid gas-dynamic problems). Also let 

L- , be the matrix L evaluated at Uw*!/, i.e., 



L tv--. 


. } 
u, / 


and 


/ 





r 

1 









( 2 ) 


( 3 ) 
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Then Harten's explicit conservative difference scheme can be written as fol- 
lows: 



VI ^ 1 y 

v:- - u- 

' _ 

AX 1 '"f+/2 

A’' \ 


(4-a) 

with 
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■ T^'Tf )'■ 
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(4-b; 

where 

A t Is the time Increment, Ax is the spatial increment, 

with 


i 










;4C 

and 

4, = 

YMJ/ ( 



(4i) 




4 ■ X •' 

’^v.i 1 j 








(4-e) 


% = 

C ^ » /TI4./: i Oj 

^ ?J-^' 




-jUa" 

“i'' 



(49) 
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For a discussion of the choice and properties of the function Q( V), the 
reader should refer to Marten's original paper. Two choices of Q that we 
have Investigated are: 


(a) Q(y) - 1/4 


(4J) 


(b) Q(.) 




J 


V > 2 f 

6 a fixed constant 


(Ak) 


The last term in equation (4b) is designed in such a way that the approxima- 
tion satisfies the guidelines listed in the previous section; that is, these 
are terms which ensure that a true physical weak solution can be obtained 
with high accuracy. 


We can view this conservative differencing as a form of the second-order 
accurate Lax-Wendrof f scheme plus an additional term. The global accuracy of 
this scheme is second-order and it is a nonlinear differencing scheme even if 
the Jacobian A is a constant matrix. In addition, Marten's method, like 
Lax-Wendrof f , has n steady-state dependence on <At. 


III. EXTENSIONS 


It is well-known that the time step of conventional explicit schemes are 
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limited by a CFL number of order one and result in long computation times for 
"stiff" problems or for steady~state calculations. For stiff problems we 
would like to retain Harten's high resolution feature but modify the scheme 
to be implicit. For steady-state calculations, we consider the following 
ways of incorporating high-resolution schemes: 

(a) Compute an approximate steady-state solution by some efficient 
scheme (explicit, implicit or a mixture of the two) and then apply an ex- 
plicit high-resolution scheme like Marten's as a post-processor. 

(b) Modify an implicit scheme in delta form [7], where the Iteration 
matrix on the left-hand side is a simple modification of the original impli- 
cit operator, while the right-hand side is the appropriate representation of 
the explicit high-resolution scheme. 

(c) Develop a fully implicit version of the high-resolution method with 
a proper linearization for the implicit operator. 

One-Dimensional Problem ; 

Preliminary numerical experiments with techniques (a) and (b) for the 
quasl-one-dlmenslonal nozzle problems showed encouraging results (see section 
V for details). However, preliminary analysis of technique (c) shows that 
Marten's scheme is highly nonlinear and special techniques need to be devel- 
oped for a fully implicit efficient Implementation; therefore, technique (c) 
will not be dlscvssed further in this paper. Mere we will describe technique 
(b) for the one-dimensional problems. 




1 
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For ilmplleityi consider Che dees of iapllclC linear one-etep tioe dif-* 
ferenclng in a noniterative delta-form [8] for the ayaten of conservative 
laws( 1) 


^ A“yu“ - 


rt, t-i 

where iU ■ U - U and should not be confused with defined in equa- 
tion (2). If 9 “ 1/2, the time differencing method is the trapezoidal for- 
mula and if d > 1, the time differencing method is backward Euler. There is 
no connection between the "6" in equation (5) with the " 9 ^" in equation 
(4d). For the spatial differencing, three-point central differencing or up- 
stream differencing is commonly used. The simplest implicit extension of 
Marten's method is to replace the right-hand side of (5) by the right-hand 
side of equation (4a). Some modification of the left-hana side of (5) can 
also be made. For comparison, this differencing of the spatial derivative 
can be viewed as three-point central differencing plus some "appropriate" 
numerical dissipation, implicit on the left-hand side, explicit on the 
right-hand side. This method will be called the "modified implicit method". 


For a steady-state calculation, the right-hand side differencing deter- 
mines the solution. We can keep the left-hand side the way it is or (to im- 
prove convergence) add a scalar dissipation term [9] to the left-hand side; 

e.g.. 


J?4X 2 AX 


( 6 ) 
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with - max • A constant version of the scalar dissl- 
pation tern of equation (6) is obtained by setting " constant for all J 
and n. This form of dissipation is sometimes called second-*order Implicit 
numerical dissipation. 

Two-Dimensional Problem! 

Although, Marten's scheme was developed for one-dimensional problems, 
the scheme can be Implemented in two spatial dimensions by applying a Strang 
type sequence of fractional step (time splitting) [10]. Again, numerical ex- 
periments for the two dimensional shock reflection problem show that the 
technique described in (<.) Is applicable. 

The extension of technique (b) in two dimensions is not straightforward. 
It is well-known that fractional step methods have the property that the 
steady-state (if one exists) depends on At even if the original 
one-dimensional version of the scheme does not depend on At. For explicit 
methods, a fractional step procedure does not introduce a serious error. 
However, if one could take a large time-step by makl.ig the method implicit, 
then the steady-state accuracy will be degraded. In general, the 
steady-state dependence on At for implicit methods can be avoided by using 
an alternating direction implicit method. Recall that Marten's 
one-dimensional scheme has an inherent dependence on At in the steady-state. 
Consequently, technique (b) will also have a steady-state dependence on At 
even if we use the alternating direction implicit method. This dilemma is 
the subject of further investigation. For the purpose of this paper, a prel- 
iminary version of technique (b) in two space dimensions is called the "modi- 
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fl«d implicit method. In the next section, we will describe the Implemente** 
tion of Harten's method and the modified implicit method for the 
two-dimensional invlscld compressible (Euler) equations of gas dynamics* 


IV. APPLICATION OF HARTEN'S METHOD TO TWO-DIMENSIONAL EQUATIONS OF GAS DYNAMICS 


In two spatial dimensions, the invlscld compressible equations of gas dy- 
namics can be written in conservative form as 


ILL j. L = 0 

St 


where 


1/ :- 


P(U)* 


\ 

I ; 


GiV) 


\v(e-¥f) 


with m •J’u and n “J*v. The primitive variables of (8) are density jo , velo- 
city components u and v, and the pressure p. The total energy per unit vo- 
lume e is related to p by the equation of state for perfect gas as 


-f -(.-■)[=- 1^1 
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where t le the ratio of specific heets. 


Let A denote Jacobian matrix 3F(U)/aU whose eigenvalues are 

( X*iX* » “ (u-c,u,u-H:,u), %diere c is the local speed of sound. The 


right eigenvectors form the matrix R,, - (R' »R^ >R|^ ) given by 




1 

1 

1 

0 

U-C 

cc 

u+c 

0 

V 

if 

\r 

1 

H -- u c 

2 


V 


H = 




Let the grid spacing be denoted by Ax and Ay such that x “ jAx and y ■ kAy. 
Using the same notation as in section II, the vector << of equation (3) for 
the x-direction (omitting the k index) is 




(da -bty'2 

(aa + bby's 


where lid 




PAGE 12 


Similarly, let B denote the Jacobian matrix 5G(U)/9U whose eigenvalues are 
/ V- \* \ 4 \ _ / 




/ 


^,v+c,v) 

. The 

right 

eigenvectors 

given 

by 



i 

i 

± 



U 

U 

LL 

1 


t --- c. 

\r 

•MC 

0 i 


M-v-C 

*■ 

2 


Ll 

> 


(13) 


The vector (X for the y-directlon (omitting j index) is 



' (cc-dd)/2 




(CC+ dcOA 




id =■ " ■'4+Kt 


where 
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As sentloned previously, the slnpllest form of 

H\i)A 05) 

Roe In [6] uses a special form of averaging Chat has Che computational advan- 
tage of perfectly resolving stationary discontinuities (for one-dimension). 
Roe's averaging takes the following form: 






'VhiC 




I (16) 




In all the computational results presented in this paper, we used Roe's aver- 
aging. 


Fract ional Step Method ; 


Marten's explicit method can be conveniently implemented in two space 
dln.3nsions by the method of fractional steps with h * ^t as follows: 




M ^ * 
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that is 


yM y>^ ^ ,n 


3;* 


tihere 






j. A< 

At 




(19) 


A 




( 20 ) 


Here, it is understood that the scalar values and the vector R'*' in the suoma- 
tions in equations (19) and (20) are values of (4c-4k) evaluated at the cor- 
responding x-coordinate by using equations (10-12) and at the y-coordinate by 
using equations (11,13,14). We use the subscripts j for x and k for y for 
the indexing of the computational mesh. Recall, for simplicity, that we om- 
itted the k index in equation (12) and omitted the J index in equation (14). 


In order to retain the original second-order time accuracy of the meth- 
od, we use a Strang type of fractional step operators, namely 


■n-‘Z 



y ^ 

ht 


( 21 ) 
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We will call equation (21) the half-step fractional step oethod. For 
steady-state calculations, the intermediate steps of equation (22) are just a 
part of the iteration procedure. Therefore if we handle the boundary data 
for the intermediate steps correctly, we can combine the adjacent 
operators into and the adjacent operators into . 

The half-step operators need only be applied at the first and last Itera- 
tions, l.e. 


it- 




We will call equation (23) the full-step fractional step method. Numerical 
experiments on the shock reflection problem show that equation (22) and (23) 
give almost identical numerical results and the amount of CPU time required 
for equation (23) is half that required for equation (22). 


Dimensions ; 


The two dimensional alternating direction implicit version of (5) is 


= - Ai: 

+6M^E)ATf = Axf 


(i^ + 


Tj"‘'= U"+AT/ 


(26) 
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- I 

i I 

■ , i 


The modified implicit method is obtslned by replacing on the 

^ right-hand side of (24) by 

j. A A 

where F and G are given by (19) and (20), and the left-hand sides of (24) and 

* 

(25) are replaced by the appropriate analogues of (6). 


^ ■' r:. 


,r 


A’* 






(27) 


V. NUMERICAL RESULTS FOR QUASI-ONE-DIMENS lONAL NOZZLE PROBLEM 

A detailed implementation of Harten's method for the one-dimensional 
inviscid compressible equations of gas dynamics can be found in Harten's ori- 
ginal paper [5]. Here, we will only describe some of our numerical results. 
For numerical experiments, we chose the quasi-one-dlmensional nozzle 

problem with two nozzle shapes (divergent and convergent-divergent). In all 
of the calculations the computational domain was 0 ^ x ^ 10. Two mesh spac- 
ing were considered for the divergent nozzle; Ax » 0.5 and Ax - 0.2. The 
I mesh soacing for the convergent-divergent nozzle was ax - 0.2. 

? Analytical Boundary Conditions ; 

^ We snecified all three primitive variables J> , u and p for the superson- 

i ■ 

ic inflow case, the two primitive variables ^ and p for the subsonic inflow 

i uase, and the variable p for the subsonic outflow case. 

if 
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f 


tf 


thmerlcal Boundary conditions : 

We use first-order space extrapolation for the outgoing characteristic 
variables [11] to obtain the nusurlcal boundary conditions for the unknown 
flow variables at the boundaries. Since Harten's scheme Is a five-point 
scheme In space, we also need the values of g> and G* on both boundaries. 
For convenience, we will use zeroth-order space extrapolation for these va- 
lues. 

The Form of 0 and : 

In all of the numerical experiments for the one -dimensional test prob- 
lems, we use equations (4j) or (4k) for the representation of the Q function 
and Roe's formula for the evaluation of Numerical experiments 

show that the two forms of the Q function have little effect on the 
steady-state numerical results of these test problems. 

Initial Conditions ; 

We use linear interpolation between the exact steady-state boundary va- 
lues as initial conditions. 

Dis cussion of Results ; 

To illustrate the shock-capturing capabilities of Harten's method and 
the modified implicit method, we compare in figures (1) and (2) the computed 
results with the first-order flux-vector splitting scheme [9], and a conven- 


k 
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tional implicit method (5) using three-point central spatial differencing 
with an added fourth-order explicit numerical dissipation [8]. From the re- 
sults, we can see a definite improvement in shock resolution by Marten's 
scheme and the modified implicit method, especially for the coarse mesh spac- 
ing. 


To demonstrate the advantage of using Marten's method as a fine tuning 
process for the conventional implicit scheme, we plot in figure (3a) a more 
difficult flow solution with the conventional implicit scheme as a 
pre-processor to obtain the steady-state solution. The undershoot and over- 
shoot near the shock are typical of the three-point central spatial scheme. 
Figure (3b) shows the improvement of the solution after we applied Marten's 
scheme as the post-processor. Figure (4) also shows the solution Improvement 
after we applied Marten's scheme as the post-processor for the solution il- 
lustrated in figure (Id). Note that it took approximately 700 steps for Mar- 
ten's method alone as opposed to 50 steps of the implicit flux-vector split- 
ting method plus 40 steps of Marten's method to converge to the steady-state 
solution. Figures (3-4) illustrate how well the method can capture the 
shock, especially for the more sensitive subsonic inflow, subsonic outflow 
case. In all of the above test problems (except for Marten's method), we 
used the backward Euler method as the time differencing. 


VI. NUMERICAL COMPUTATION OF TWO-DIMENSIONAL SHOCK REFLECTION PROBLEM 


Various first-order upstream conservative differencing methods have been 
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developed for systems of one*^ilmenslonal hyperbolic conservation laws 
[6,9,12-14]. When these methods are applied to the one-dlmenslonal Invlscld 
equations of gas dynamics, most of the methods capture shocks very well. 
However, a direct application of these methods to problems In two-dlmenslons 
does not necessarily produce the same type of high resolution near shocks un- 
less the shocks are aligned with one of the computational coordinates. This 
Is due to the fact that first-order methods are generally not adequate to 
solve a complicated flow field accurately except on a very fine mesh. On the 
other hand higher-order extensions of these first-order upstream methods are 
usually rather complicated to use (see for example, van Leer [15]). Harten's 
method Is less complicated chan other recently proposed second-order 
hlgh-resolutlon n^thods; however. It was developed for one-dlmenslonal prob- 
lems. 


In order to see how well Harten's method and the modified implicit meth- 
od can capture shocks for two-dimensional Invlscld gas-dynamic problems, we 
consider a simple invlscld flow field developed by a shock wave reflecting 
from a rigid surface (fig. 5). The steady-state solution can be calculated 
exactly and thus can aid us in evaluating Che quality of Che numerical meth- 
od. Figure 5 shows the Indexing of the computational mesh. 

Analytical Boundary Conditions; 

The boundary conditions are given as follows: (a) supersonic Inflow at 
j»l, k«l,...,K which allows the values to be fixed at freesCream 
conditions; (b) prescribed fixed values of U . at k ■ K, j «• which 
produce Che desired shock strength and shock angle; (c) supersonic outflow 
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at j >■ J, k a 1 K; (d) a rigid flat surface at k - 1, j ■ which 

can be shown to be properly represented by the condition v ■■ 0, with the ad- 
ditional condition dp/»y - 0 at k « 1 from the normal y-momentum equation. 

Initial Conditions ; 

Initially, the entire flow field is set equal to the freestream super- 
sonic Inflow values plus the analytical boundary conditions as described 
above . 


Nu merical Boundary Conditions fill ; 


The supersonic outflow values U_ , k = 1,...,K-1 are obtained by zer- 
oth-order extrapolations, i.e.. 




k = 1 K-1 


(28) 


The values of ji. , on the rigid surface with j 
ined by zeroth-order extrapolations, i.e. 


are also obta- 








(29) 


Three different methods were used in approximating ^p/dy 


0 to get e » 
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(a) normal derivative of e equal to aero (f Irat-order) 




(30a) 


(b) normal derivative of e equal to zero (aecoiid-order) 


e. - (4e , - 6 a » )/3 


(30b) 


(c) normal derivative of p equal to zero (second-order) 




Equation (30a) together with equation (29) Is an approximation to 
, ■ P'« ; l.e., a first-order approximation for the normal derivative of 

p equal to zero. Equation (30b) together with equation (29) is an approxima- 
tion to equation (30c). We use equations (30a) or (30b) for the implicit 
methods mainly because of their ease of application with implicit numerical 
boundary conditions. From the numerical experimental we found that equation 
(30b) or (30d) produce better numerical solutions near the wall than equation 
(30a). 


Since Harten's scheme Is a S-polnt scheme (In each spatial direction), 


we also need the values of g , 
^ I for j * 1|...|«7 or h * li...,E. 




f 

f 

1 


’ S)K ’ %)1 ’ 

For convenience, we will set 

(Si) 

%K = e,- g.i 


7 
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The Form of 0 and 


In all of the numerical experiments for the shock reflection problem, we 
use equation (4j) for the representation of the Q function and equation (16), 
(Roe's formula) for the evaluation of 

Boundary Data for the Intermediate Steps t 

When fractional steps are used for the two-dimensional problems, there 
are Intermediate boundary data which are not required for non-fractlonal step 
and alternating direction Implicit methods. For example. In order to advance 
from equation (17) (the operator) to equation (18) (the <?(^ operator), we 
need values of for j - 1,...,J. Since we are solving for steady-state 


solutions, we set v^ 




0 “ Vj' £ St the rigid surface, and j ■ 1 J 


equal to the precribed boundary values. For the remaining Intermediate boun- 

•A A A 

dary data P. , m > . and e>^ , we use equation (17) to solve for the values. 

The reason is that we know the values of the entire flow field at level n. 
Therefore we can advance to the next step by equation (18), and update the 
numerical boundary conditions by equations (28-30). We can repeat the pro- 
cess by using the sequence Indicated In equation (22). 

For the full-step fractional step method, every step Is an Intermediate 
step. The process is as follows: 


step 1: the operator 


r"* 




2 AX 


( 32 ) 
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rht value of F evaluated at the Initial value of 


%>• 


y-fi 

Step 2: the operator 


■= 




/V ■# \ 


(33) 


step 3: the o(y operator 


I 


V, - U', 
Ik 


, / A A ^ ' 

ft, _ -p. 


(34) 


Thus, we can start out by setting the intermediate values of v* 




Ji 




j * !,•••, J 


the sane way as 


yYi 

the half-step method at step 1 ( operator). The reason is that we know 
the values of the entire flow field by the initial conditions and the analyt- 
ical boundary conditions. Therefore we can advance to the next step by equa- 

•• 

tion (33). After the values of Uv . , k - 1,...,K and j * 2,...,J are calcu- 

lated, we use the precrlbed values at the inflow as U .and equation (28) as 
U ^ In equation (33) to solve for U. . and U*. . At this stage step 2 

TyK. 

( Operator) is finished and we can apply similar process to step 3 

( operator). We can repeat the process by using the sequence shown in 


equation (23). 


Modified Implicit Method; 
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Recall we use of the term "modified Implicit method" to denote the al- 
gorithm obtained by applying technique (b) of section III to equation (7) 
(l.e., by the alternating direction method as described at the end of section 
IV), A fully Implicit form of the Marten's method results In a more compli- 
cated solution Inversion procedure. Work Is underway to develop a proper li- 
nearization procedure. Here we only study the applicability of the less com- 
plicated iteration matrix as the implicit operator. 

Other Numerical Methods; 

We choose four methods for comparison with Marten's method and the modi- 
fied Implicit method: (a) a first-order explicit method of Roe [6], (b) the 
first- and second-order flux-vector splitting method [9], (c) the 
explicit-implicit MacCormack's method (finite volume) [16], and (d) a conven- 
tional Implicit method with three-point central differencing in space and a 
fourth-order explicit numerical dissipation term [8]. In the calculations 
with the implicit methods (except MacCormack's method), we use the backward 
Euler for the time differencing. 

Discussion o f R esults t 

The main purpose of the two-dimensional numerical experiments Is to eva- 
luate the shock resolution of Marten's method and the modified Implicit meth- 
od. Convergence rate and computational efficiency will oe investigated and 
reported elsewhere. We expect that techniques (b) and (c) proposed in sec- 
tion III will speed up convergence for subsonic (dominated) steady-state 
problems. 
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In all of the numerical experiments, the Incident shock angle Y* was 
o 

29 and the freestream Mach number was 2.9. The cunputatlonel domain was 

0 ^ X 4.1, and 0 ^ y ^ 1, and the grid size was 61X21. The appropriate 

analytical boundary conditions are applied along the boundaries cf the doma- 
in. Only pressure contour and pressure coefficients will be Illustrated. 

2 

Here, the pressure coefficient is defined as c* ■ 

r 

p the freestream pressure. The exact minimum pressure corresponding to 

•O 

'\f'm 29* and “ 2.9 Is 0.714286 and the exact maximum pressure Is 2.93398 

(see Fig. 6). Forty-one pressure contour levels between the values of 0 and 
4 with uniform Increment 0.1 were used for the contour plots. The pressure 
coefficient was evaluated at y - 0.5 for 0 < x < 4.1. 

The exact pressure solution Is shown In figure 6. Pressure contours and 
the pressure coefficient evaluated at y “ 0.5 are shown in figure 7 for seven 
different methods. Both the first-order explicit method of Roe and the im- 
plicit first-order flux-vector splitting method yield smeared discontinuities 
as shown In figures (7a) and (7b). The second-order (in space) flux-vector 
splitting method, the second-order (In space) conventlonax Implicit method, 
and MacCormack's explicit-implicit method are shown in figures (7c)-(7e). 
(The results for MacCormack's method were provided by W. Kordulla, National 
Research Council resident research associate at NASA Ames Research Center.) 
Harten's method with full-step fractional steps and the modified implicit 
method are shown In figures (7f) and (7g). These methods show a definite Im- 
provement in shock resolution. There are practically no oscillations over 
the entire flow field. Figure 8 shows the pressure contours of Harten's 
method by using the halT-step fractional step method. The solution Is nearly 
Identical to figure (7f) but required half the number of iterations (l.e.. 


- i) with 
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half of the step size, double the CPU tiae). Figure (9a) shovs the results 
after 30 steps of the conventional implicit method. Figure (9b) shows the 
improvement after we applied 200 steps of Hartea's method as post processor. 
All of the above methods required 15''-600 iterations to converge with CFL 
numbers ranging from 0.5 to 0.8 (except figure (9a)). 

At the present stage of development, the prllminary version of the modi'- 
fled implicit method for two-dimensional problems has several deficiencies. 
In particular, the Implicit operators (left-hand sides of (24) and (25)) are 
not proper linearizations of the right-hand side. In fact, the lineariza- 
tions are crude at best, and one would expect chat this would have a serious 
effect on the stability of the method. In addition. Marten's method, like 
Lax-Wendroff , has a steady-state dependence on ^t. For explicit methods, 
this does not introduce a serious error. However, if one could take a large 
time-step by making the method implicit, then the steady-state accuracy is 
expected to be degraded. Figure 10 shows the pressure contour for the shock 
reflection problem at a CFL number of 1.2. Even at this modest CFL value, 
the accuracy is degraded as compared to the same problem run at a CFL number 
of 0.5 shown in figure 7g. 


VII. CONCLUSION 

The application of Marten's explicit method to quasi-one-dimensional 
nozzle problems and to a -dimensional shock reflection problem resulted in 
high shock-resolution steady-state numerical solutions. By combining the ad- 
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Jacent half-step operators Into full-step operators of the fractional step 
nethod, one can cut the computation time in half. Applications of the 
post-processor method and the modified implicit method for steady-state cal- 
culations show encouraging results for one-dimensional problems; however, 
testing in two dimensions is not complete and further investigation is needed 
for efficient implementation of the implicit method. 
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(a) Hartea's aethod (explicit). 


(b) Hodlfied implicit aethod. 




(c) Conventional implicit aethod. 


(d) First-order lapllcit flux-vector 
splitting aethod. 


Fig. 1 Density distribution: supersonic inflow, subsonic outflow. 

Nozzle shape: " ; 20 spatial Intervals. 
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(a) Harten's oethod (explicit). 



(b) Modified Inpliclt nethod. 



(c) Conventional implicit method. 


(d) First-order implicit flux-vector 
splitting method. 


Fig. 2 Density distribution; supersonic inflow, subsonic outflow. 


Nozzle shape: 


9 


50 spatial intervals 
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(a) Conventional implicit method. (b) Batten's method as post-processor. 


Pig. 3 


Density distribution: subsonic inflow, subsonic outflow. 


Nozzle shape: 



50 spatial intervals. 






Fig. 4 Density distribution: supersonic inflow, subsonic outflow. 

Batten's method as post-processor (use the result of fig. 
(Id) as initial conditions). 







Pig* 5 Indexing of conputational mesh 
for shock reflection problem. 
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(«) 30 steps of the conventional (b) 200 steps of Harten's nethod 

Implicit aethod. ^be result of fig. (9s) 

as initial conditions). 

9 Pressure Contours* Harten's aethod as post**processor. 




Fig. 10 Pressure Contours: Modified 

iaplicit aethod at CFL “ 1*2 
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